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ABSTRACT 

We propose a new, physically motivated fitting function for density PDFs in turbulent, ideal gas. 
Although it is generally known that when gas is isothermal, the PDF is approximately lognormal in the 
core, high-resolution simulations show large deviations from exact log-normality. The proposed function 
provides an extraordinarily accurate description of the density PDFs in simulations with Mach numbers 
~ 0.1 — 15 and dispersion in log (p) from ~ 0.01 — 4dex. Compared to a lognormal or lognormal-skew- 
kurtosis model, the fits are improved by orders of magnitude in the wings of the distribution (with fewer 
free parameters). This is true in simulations using a variety of distinct numerical methods, including or 
excluding magnetic fields. Deviations from lognormality are represented by a parameter T that appears to 
increase systematically with the compressive Mach number of the simulations. The proposed distribution 
can be derived from intermittent cascade models of the longitudinal (compressive) velocity differences, 
which should be directly related to density fluctuations, and we also provide a simple interpretation of 
the density PDF as the product of a continuous-time relaxation process. As such this parameter T is 
consistent with the same single parameter needed to explain the (intermittent) velocity structure functions; 
its behavior is consistent with turbulence becoming more intermittent as it becomes more dominated by 
strong shocks. It provides a new and unique probe of the role of intermittency in the density (not just 
velocity) structure of turbulence. We show that this naturally explains some apparent contradictory results 
in the literature (for example, in the dispersion-Mach number relation) based on use of different moments 
of the density PDF, as well as differences based on whether volume-weighted or mass-weighted quantities 
are measured. We show how these are fundamentally related to the fact that mass conservation requires 
violations of log-normal statistics. 

Key words: hydrodynamics — instabilities — turbulence — star formation: general — galaxies: forma- 
tion — galaxies: evolution — galaxies: active — cosmology: theory 



1 INTRODUCTION 

It is well-established that idealized isothermal turbulence (i.e. 
driven turbulence in the absence of external forces) drives the den- 
sity probability distribution function (PDF) to an approximately 
log-normal shape (a Gaussian in lnp), with a dispersion that in- 
creases weakly with the Mach number (see e.g. Vazquez-Semadeni 



[T994l|Padoan et al.p997l |Scalo et al.|1998||Ostriker et al.|1999| 
Klessen 2000). This is easily understood: the Navier-Stokes equa- 
tions for inviscid, isothermal gas can be written only in terms of 
lnp (i.e. changes due to the velocity field are multiplicative), so 
if the density distribution is the product of a sufficient number of 
random, uncorrected and unbiased multiplicative perturbations it 
should converge via the central limit theorem to a lognormal ( |Pas-| 
|sot & Vazquez-Semadeni|[T998l |Nordlund & Padoanp999) . This 
assumption has become widespread, and forms the basis for inter- 
pretation of many observations (e.g. Ossenkopf & Mac Low|2002| 
|Ridge et al. 2006] [Wong et al.||2008| , as well as a wide range of 
theoretical models for ISM structure and star formation ( Elmegreen 
|2002l|Pado"an & Nordlund|2002| [Krumholz & McKee|2005HHerTi 
|nebelle & Chabrier||2008"[ |Hopkins|[20TT| |2012a|c|d| Federrath & 
Klessen in preparation). 

However, a growing number of higher-resolution simula- 
tions show deviations from exactly log-normal behavior in driven 
isothermal turbulence, which can become very large (orders-of- 
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magnitude) in the tails of the distributions (see e.g. Federrath et al. 
[20081 |20T0l|Federrath & Klessen|20i2[|Schmidt et al.|2009||Price 
& Federrath 2010; Konstandin et al. 2012). This appears to be es- 
pecially severe when the density PDFs are broad (high- variance). 
And there are well-known, highly statistically significant inconsis- 
tencies between the properties inferred from the volume-weighted 
versus mass-weighted density PDFs, or between different moments 
of the density PDFs, which should not appear if they were truly 
log-normal (see above and Lemaster & Stone 2009 ; Burkhart et al. 
[200^ [Federrath et al.|2010[|Price et al.|2011[|Molina et al.|2012] K 

In fact, mass conservation fundamentally invalidates the as- 
sumptions that lead to the log-normal prediction. First, consider the 
density distribution sampled within a large box, measured as the 
volume-average density p, = m,/vo in independent, equal, small 
sub-volumes vo, and assume it is exactly lognormal. The density 
distribution smoothed on a larger scale vi 2> vo is then necessar- 
ily the sum over the sub- volumes vo. , inside vi , p, = m,) /vi = 
(vo/vi) Pu i- e - it is m e linear sum of lognormally distributed 
variables. But the convolution of linear values of lognormal vari- 
ates is not lognormal (although it can be approximately so; see 
|Fenton|[l960} . So mass conservation implies the density distribu- 
tion cannot be lognormal on more than one scale. Second, conver- 
gence to a lognormal PDF via the central limit theorem requires 
fluctuations be uncorrelated; but this cannot be true in detail. Con- 
sider again a small but finite volume vo within some larger volume 
containing total mass M; there is nothing stopping there from be- 
ing arbitrarily small mass within vo (i.e. fluctuations can extend 
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Figure 1. Examples of the proposed (volume-weighted) log-density PDF 
(Eq.[5j, as a function of parameters T (increasing with deviations from log- 
normal) and Si„p, v (volume-weighted variance in lnp). Top: Varying T at 
fixed variance. T = is lognormal (the black solid line); the skew increases 
systematically with T. Bottom: Changing variance at fixed T (broader PDFs 
at larger S) np y). A simulated box averaged over progressively larger scales 
should exhibit approximately fixed T with decreasing Si„ p . v , varying the 
PDF shape as shown. 



to infinitely large in negative lnp), but mass conservation sets an 
upper limit to the maximum volume-averaged density within the 
volume (lnp < ln(M/vo); the local density on smaller scales is 
free to diverge). So especially for large fluctuations, which may 
be important in the tails of the distribution, mass conservation im- 
plies a biased/asymmetric distribution in lnp. There are many other 
reasons the PDF may be non-lognormal; for example, it is well- 
known that compressive (longitudinal) velocity fluctuations are not 
exactly Gaussian-distributed, owing to intermittency (discrete phe- 
nomena such as shocks or waves producing extended tails in the 
PDF), and there are (small) residual correlations between the local 
Mach number and density which imply correlated fluctuations (see 
Passot & Vazquez-Semadeni 1998; Kritsuk et al. 2007; Federrath 
et al.|2010) . " 

This is important for several reasons. Many of the physically 
interesting behaviors in turbulence, and models described above, 
rely critically on the behavior in the extremes of the distribution 
(where deviations might be very large). Also, the assumption that 
PDFs are lognormal, and hence that their moments are related by 
certain simple analytic scalings, has led to widespread use of "prox- 
ies" for other quantities of interest (for example, it is common 
to infer, rather than measure, the dispersion in ln(p), by actually 
measuring the mean value, and assuming such a linking relation) - 
these could be seriously in error. This could severely bias estimates 
of the dispersion-Mach number relation in different studies. Non- 
normal behavior in density PDFs may provide a valuable additional 
probe of intermittency in compressible turbulence, which remains 
poorly understood. And the differences from lognormal behavior, 
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Figure 2. Fits to simulated isothermal turbulent boxes. Top: Volume 
weighted PDFs from simulations (points) in Table ^ spanning a range 
in Mach numbers and variances, driving schemes, numerical methods, 
and magnetic field strengths. In the same color, lines show the best-fit 
with the proposed Eq. [5] Bottom: One example (the |Federrath et al.|2010| 
compressively-forced M = 5.6 model), comparing the best-fit with Eq.[5]to 
a standard lognormal fit, a lognormal fit with independently freed normal- 
ization and mean (i.e. not constrained to be properly normalized or conserve 
mass), and a lognormal-skew-kurtosis model (fourth order log-log polyno- 
mial). The proposed model fits extremely well, whereas even higher-order 
lognormal models (with more free parameters) fare very poorly (especially 
in the tails of the PDF). 



especially in how the volume-weighted and mass-weighted PDFs 
relate, are critical to relate results from different simulations to one 
another and to understand convergence in different numerical meth- 
ods. 

In this paper, we propose a new functional form for the den- 
sity distribution in isothermal, driven turbulence. We show that this 
provides an excellent fit to the density PDF in simulations (signif- 
icantly improved over a lognormal model), over a wide dynamic 
range. We also show that this proposed PDF resolves the incon- 
sistencies with mass conservation discussed above, and can be de- 
rived from models of intermittency in turbulent velocity fluctua- 
tions, linking the non-normal features in density PDFs to funda- 
mental properties of compressible turbulence. 



2 MODELS FOR THE DENSITY PDF 

Consider the density PDF. We define first the volume-weighted PDF 
of log-density, fV(lnp), i.e. the probability of a given unit volume 
V having an averaged density p = M(V)/V, per unit In p. This must 



© 0000 RAS, MNRAS 000, 000-000 



Intermittent Density PDFs 3 



1.0 



T~0.1 s lnp , M 

T~ 0.25 ln(1 +0.25 Sf np , M ) 




0.01 0.1 1 10 

Compressive Mach Number af c 

Figure 3. Top: Fitted parameter T representing the deviations from log- 
normal, as a function of mass-weighted variance in In p, for each simu- 
lation in Table [T] T appears to increase with Si np .M, m simple fashion. 
Bottom: Same, but versus the mean compressive component of the driv- 
ing Mach number Ai c (Mc ~ M. for pure compressive driving, pb M/3 
for solenoidal driving). The dependence is even more pronounced, with T 
rising sharply at Jv[ c > 1, approximately as T ~ 0.05 M c . 



(1) 



obey two constraints: 

/oo 
Py(lnp)dlnp = 1 
-OO 

/OO 
p/V(lnp)dlnp = p = 1 
- OO 



the first condition is just that the PDF be properly normalized, the 
second is mass conservation (integrated over all volumes V, the to- 
tal mass M is conserved: J pV M Pv dlnp = M = po Vi 0l ). Through- 
out this paper, we use dimensionless units such that the volume 
averaged mean density po = 1 . 

We can also define the mass-weighted density PDF, Pu(lnp). 
Since p = dM/dV, it follows that f M (lnp) = pfV(lnp) (see e.g.O 
|et al.|2003| >. 

2.1 Lognormal PDFs 

In the simple log-normal case, we have: 



/V LN (mp) = 



1 



\/27r5i„p.\ 



: exp 



(lnp + 5in P ,v/2) 2 

2Sl n p,v 



(2) 



where 5i n p.v is the volume- weighted variance in lnp (Sinp.v = 
J(lnp— (lnp)) 2 Pv(lnp)dlnp). The normalization and mean are 
fixed by Eq.[T| In this case, the various volume- weighted moments 




lnp, M 



Figure 4. Moments of the density distribution. The simulations from Ta- 
ble[T]are shown, with the prediction for a lognormal model (Eqs. |3|4) , and 
our proposed Eq.|5](see Eqs. |6|7( f or fixed T = 0.05, T = 0.5 and the best 
fit T = T(S) oc M c from Fig. |3| Top: Mass-weighted dispersion in lnp 
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) versus volume-weighted dispersion in lnp (= 



: sV 2 ,,). These 

In p. M' 



are equal in a lognormal model. The simulations do not follow this relation, 
but trace the mean T = T{S) track between T = 0.05 and T = 0.5. Middle 
Volume-weighted dispersion in the linear p, versus volume-weighted dis- 
persion in lnp. The simulations lie well below the lognormal prediction. 
Deviations from lognormal are significant at all S here (though difficult to 
see at smaller values), but agree well with the proposed model. Bottom: 
Volume-weighted dispersion in the linear p, versus mass-weighted disper- 
sion in lnp. Although the data still favor the non-lognormal model, it is 
worth noting that deviations from the lognormal case are minimized when 
comparing these specific moments. 



of p and lnp are easily calculated: 



(lnp)y N = — Sinp.v/2 



(3) 



Jp,V 



: exp (Si 



In p, V ) 
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The mass-weighted PDF PM(lnp) is also lognormal, with moments 



cLN _ c 
>Jliip,M — Jlnp, V 

iLN 



(lnp) M = +S{ np ,M/2 

(p)Jf = exp(S] np , M ) 



(4) 



Jp.M 



= exp(35i„ p .M) -exp(25i„ 



p.Mj 



The quantity T clearly represents a degree of non- 
lognormality, and breaks the similarity otherwise seen between e.g. 
the volume-weighted and mass-weighted distributions. Note that as 
T — > 0, these relations do simplify to those of the lognormal distri- 
bution; Eq.[5]becomes identical to Eq.|2] Fig.[T|illustrates the shape 
of this PDF for different T and 5i np .v. 



2.2 An Improved Model 

Motivated by the considerations in §|4]& Appendix|A] consider the 
following simple model for the volumetric density PDF: 



(5) 



Pv (lap) dlnp = h (2 V A u) exp [- (A + u)] \ - du 



m=0 
A 



A e 



(m-1) 



■ d;/ 



A = 



i + r 

Slnp.V 

2T 2 



hip 



(«>o) 



where /i(jc) is the modified Bessel function of the first kindQand 
Pv — for u < 0. It is straightforward to verify that this obeys both 
constraints in Eq. [T] Given some variance S\ np ,v (or equivalently, 
parameter A), then, this is a one-parameter model in the value of 
the adjustable parameter T. We will discuss the meaning of this 
parameter below. 

From this, we can immediately derive the following useful re- 
lations for the moments: 



S{np,V - 

Sp.v - 

(lnp) v ■■ 

(P)v- 

S\n p.M : 



2XT 
■ exp 
; exp 



(6) 



2Ar 

1 + 3 T + 2 T 2 

Sinp.V 



1+3T + 2T 2 

xt 2 S] np ,v 



l+T 



2 XT 

(i+fy 



i 
i 

(i + r) 



= Si np , v (l+T) 



(V) 



Sp, M = exp 



6 XT 2 



1+4T + 3T 2 



exp 



4 XT 



I 35i„ p .m (l + T 
:exp(^ — 



(In p) M ■- 

(p)m ■- 



XT 1 



4T + 3T 2 



1 + 3T + 2T 2 

/ 25i„ p .M(l + r) 3 
6XP l 1+3T + 2T 2 



(l+T) 



(l + T) 2 2 

;exp (iww) =exp ( 



S]np,M (1 + T) 3 

1 + 3T + 2T 2 



3 RESULTS 

Fig. [2] compares our proposed Eq. [5] to an ensemble of simula- 
tion results compiled form the literature. Each simulation is an 
idealized driven "turbulent box," but we intentionally consider a 
wide range of results using different turbulent forcing schemes 
(solenoidal vs. compressive), strength/Mach numbers (from M ~ 
0.1 — 15), numerical methods (fixed grid, adaptive-mesh refine- 
ment, and smooth-particle hydrodynamics), and resolution. The 
properties of each simulation are summarized in Table [T] For our 
purposes this allows us to examine the PDF over a wide dynamic 
range. In each case we take the volume-weighted PDF Py (or con- 
vert P M presented by the authors to P v ). 

We fit each simulated PDF taking S\ np ,v as fixed to the true 
variance in each distribution]^ so that there is only one free param- 
eter T. In every case, our proposed function fits the data excep- 
tionally well; there is no case with % 2 l v 2> l|^]ln contrast, in a 
majority of cases (and every case with S\ np , v ^1) the lognormal 
model is inaccurate (x 2 / u 2> !)• 

This is illustrated for one typical PDF, with a large variance 
Sinp.v ~ 9.2. We show the simulation data and our best-fit with 
T w 0.4. The RMS deviation <5 rms of the simulation data about the 
best-fit is < 0.05 dex. In contrast, the lognormal model with the 
correct variance S\ np .v has « 0.77 dex. If we free all param- 
eters in the lognormal model, so that it now has three degrees of 
freedom (as opposed to one in our proposed model), the best-fit 
Sims = 0.51 dex (moreover, this model produces a PDF which is 
non-unitary by a factor of ~ 5, fails at mass conservation by a factor 
of w 80, and under-predicts the correct variance by ~ 1.2 dex). We 
could simply consider higher-order polynomials in log-log space: 
for example, fitting a lognormal-skew-kurtosis model (i.e. adding 
third and fourth-order terms in the exponential as proposed in |Fed-| 
|errath et al.| pOTo) ); but even with these added degrees of free- 
dom we find the best-fit has 5 rras = 0.4 dex. In fact, provided that 
the necessary physical constraints (Eq. [T]) are obeyed, we require 
an eighth-oidei log-log polynomial to recover the same accuracy 
in <5rms as our proposed fit. This not only introduces seven addi- 
tional free parameters, but it also is extremely unstable (the results 
can change by orders of magnitude if we add/remove points in the 
wings of the distribution). 

Fig. [3] shows the magnitude of T (the strength of deviations 
from lognormal statistics) fit to each simulation, as a function of 



I\ (x) satisfies 



(x/2) 



2m+l 



„ m! (m+ 1)! 

2 d 2 h 9/j . 

x — - + x (x +l)h = 

dx 2 dx 

A useful relation for large x > 100 is 

h(x) -* (2nxy' /2 exp(.t) 



2 We define the "true" variance as the variance measured directly from the 
simulation density PDF data points (either mass or volume-weighted as ap- 
propriate), without reference to any fitting function. 

3 The fits are also remarkably stable. If, for example, we also free Si nPi y, 
we recover nearly identical answers (well within the 1 a error bars), albeit 
now with two free parameters. We can in principle also free the normaliza- 
tion and value of po (though this would violate the constraint Eq.[TJ and 
recover the same answer within ~ 1%. We have also experimented with 
randomly removing 1, 10, 30, 50 and 75% of the data points in every PDF, 
and find in every case that the results remain consistent with their \a uncer- 
tainties. 
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variance. There is a clear trend, where the deviations are small (al- 
though we emphasize in many cases still extremely statistically sig- 
nificant) at small S then rise with increasing 5. We discuss possible 
reasons for this below, but for now simply note it can be param- 
eterized by a simple linear fit T ~ 0-lS\n P ,M, or, more accurately, 
as proportional to the compressive component of the Mach number 
M c (see Table[T}, T ~ 0.05 M c - 

Fig.[4]shows how the proposed model can explain many devi- 
ations from lognormal statistics seen in the moments of these dis- 
tributions. First, we compare the standard deviations (= S 1 / 2 ) of 
lnp in both mass and volume- weighted measures. If the PDFs were 
lognormal, S\ up .m = S\„ p ,v', but we see this clearly is not the case 
(the lognormal prediction is ruled out as the null hypothesis here 
at > 6ct confidence). Simulations always find 5i np ,M < Sin P ,v, with 
the differences being especially large at high variance. We compare 
the predicted relations for T = 0.05 and T — 0.5, which bracket the 
data at low and high variance, respectively. We also compare the 
track for the mean linear relation T oc M c , which runs through the 
locus of points. 

We next compare the (volume-weighted) dispersion in the lin- 

1 12 

ear density (S v ) to that in the logarithmic density. A lognormal 
model predicts S p ,v = exp(5i np ,i/) — 1; this is again clearly ruled 
out (> 8 a confidence). Again simulations lie between the predicted 
relations for T = 0.05 and T = 0.5, with the mean T — S track 
predicting a curve that goes through essentially all points. Interest- 
ingly, if we compare S' v to the mass-weighted dispersion in the 
log density Si„p,M, it turns out that the deviations from the lognor- 
mal prediction with T are minimized. The data still favor a model 
with T > 0, but the predicted tracks differ by a small amount. This 
explains why certain properties appear to agree better with lognor- 
mal models when taken as mass-weighted (discussed below; see 
also|Konstandin et al.|2012) . 

4 A PHYSICAL MODEL FOR THE DENSITY PDF 
4.1 Quantized Log-Poisson Cascades 

We now consider a physical motivation for Eq.[5] Perhaps the most 
popular model for the statistics of intermittent turbulence is the 
| She & Lev eque (1994) (hereafter SL) cascade model. In strictly 
self-similar (Kolmogorov) turbulence, the velocity structure func- 
tionsS p (/f) = {6v(R) p ) = (|v(x) -v(x+R)p) scale as power-laws 
S P (R) oc R^ p with f p = p/3. The SL model proposed an alternative 
phenomenological scaling: 



(8) 



with the original choices /3 = 7 = 2/3 (giving Q p — p/9 + 2[l — 

(2/3)"/ 3 ]). 

|She & Waynure] {1995} and |Dubrulle| < fT994) > showed that 
the scaling in Eq. [8] was the exact result of a class of quan- 
tized log-Poisson statistics. They assume extended self-similarity, 
i.e. Sv(R)" oc R" /3 e£ /3 (where e R is the dissipation term between 
scales), and a general hierarchical symmetry between scales - such 
that the statistics of can be encapsulated in a term ttr such that 
eR ~ ttr (where describes the scaling of the average proper- 
ties of the most extreme/singular objects, e.g. shocks). Under these 
conditions, SL-like scalings are equivalent to the statement that tvr 
obeys log-Poisson statistics of the form 



P{tv r ) dn R = dY P\W (m) Gr(Y, m), Y . 
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Figure 5. Structure functions and their (potential) relation to non-lognormal 
density PDFs. Top: Structure function scalings C, p (from (Sv(R) p ) oc R^p) 
versus order p. We compare the Kolmogorov prediction £ p = p/3 (no in- 
termittency), to the She & Leveque 1 1994 1 "quantized log-Poisson" model 
for intermittency (Eq. [8] with the labeled p, 7) and the |Boldyrev| 
revision of this model for super-sonic turbulence. We also compare the pre- 
diction from Castaing 1 1996) for the "thermodynamic" model (a continuous 
compound log-Poisson-exponential distribution), which predicts a PDF for 
longitudinal velocity moments that motivates our proposed density PDF, 
with a similar parameter T v that applies to the intermittency in the veloc- 
ity field (Eq. |14| for T v chosen to match the same assumptions as the other 
intermittency cases). The two produce nearly identical structure functions. 
Data from |Pan et aT|(2009) from one of the simulations we study (see Ta- 
ble HI is shown. Middle: Structure functions £ p from the simulations in 
Fig. |2] (points), where available, and best-fit T„ intermittency model. Bot- 
tom: Comparison of T fit to the lnp distribution (Eq.|5J and T v fit to inter- 
mittency in the structure functions (Eq. 1 14) . The two increase together at 
large values, though there may be a "floor" in T v that is difficult to measure 
from the density PDF. 
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Figure 6. Dispersion-Mach number relation, with the simulations from Ta- 
ble [T] (only including those not magnetically-dominated). Top: Volume- 
weighted dispersion in lnp, versus compressive Mach number Jv[ c (see 
Fig. |5J. Though a relation is clear, it does not follow the usual 5 = 
ln(l + A^l) curve (closer to Li2(A^^), the dilogarithn0), and has large 
scatter at high-.Mf The scatter becomes much larger if we just consider M 
instead of M c . Bottom: Same, but using the mass-weighted dispersion in 
In p. Here the data agree very well with the predicted relation, with smaller 



scatter. This is related to the fact that the S\ nPt M - 
deviations from lognormal (i.e. allows for S p y • 
rections from non-lognormal terms). 



with 



S p y relation minimizes 
' Ail with minimal cor- 



PQnp) 



with /3 p = /3 1/3 (see 



1 \ - , . / In p 
— ^P MR)H G R [ — , 



(11) 



i Ho"pkins||2012bt ). 
In [She & Waymire| ( |1995| >, they assume the simplest possi- 
ble case, that Gr(Y, m) — > 5(Y — m) (a Dirac 5-function), which 
corresponds to all driving "events" (shocks, rarefactions, vortices) 
which produce velocity or density changes being strictly quantized 
(i.e. producing an exactly identical multiplicative effect in every 
instance). This would predict /'(lnp) — > P\(m = lnp/ln/? p ). In 
such a quantized log-Poisson process, the density on a small scale 
R would be related to the mean density on the box scale L by 
p(R) = ft™ (R/L) 1 po where in is Poisson-distributed with integer 
values (per P\(m)). While this can (for the appropriate choice of 
P and 7) give structure functions in good agreement with experi- 
ment, it is too simple a model for the density PDF because it pre- 
dicts In (p) itself is quantized in units of A lnp = ln/3 p , rather than 
being a continuous variable. 

4.2 Continuous Models for Intermittency 

A more general statistical cascade model which produces continu- 
ous distributions is the "thermodynamic" model proposed in |Cas-| 
taing (1996). The initial assumptions there are quite different from 
SL (for example, the scaling of the "most singular structure" does 
not need to be assumed). However, as shown in |He et ak] ( |1998| >, 
this is just a more general consistent limit of the same hierarchy 
(different choice of Gr). This model predicts directly the scaling of 
the longitudinal velocity increments ((\8v\ p ) oc a''), rather than the 
indirectly related quantities -kr or cr, and models their scaling as: 



P(v)dv = £ 



\ m e- 



m\ (m — 1 ) ! 



dv 



(12) 



where v = T v ~ 1 1 In (ov/<tl) | with T v a constant. 

We see immediately that this has the form of our proposed 
scaling for the form of A- (lnp) (with T v playing an identical role 
to T), with the assumption of Eq.[TT|that we can make a one-to- 
one mapping between the statistics of the magnitudes in the log- 
velocity field and the log-density field - i.e. that something like a 
dispersion-Mach number relation is universal across scales. 

But it is also straightforward to see that this is just the general 
form of the log-Poisson statistics with 



P\(m) 



exp(— A) 



(10) 



and Gr is any well-defined, infinitely divisible probability distri- 
bution function (physically depending on the driving and character 
of "structures" in the turbulence). This describes a general class of 
random multiplicative processes that obey certain basic symmetry 
properties. 

Since 8v(R) oc ej/ 3 oc nj 3 , In (5v/(8v)) = ln7r 1/3 = (l/3)ln7r 
is a linear transformation and should obey the same statistics as 
(1/3) ln7r R (see e.g. |She & Waymire|[l995||Dubrulle|1994| ). And 
since the basic derivations of the relation between the density and 
velocity power spectra are based on the (linear) density field being 
the product of compressive modes in the velocity field (e.g. the 
dispersion-Mach number relation), Sp = p/pa — 1 obeys the same 
statistics as 8v(R) (for small "steps" in scale). Thus the statistics of 
lnp should have the same form as (1/3) In tvr ; and from Eqn.|9] the 
statistics of (1/3) ln7rj; are identical to the statistics of Iutvr for a 



G R (Y\m)dY 



exp (— v) dv 



(w-1)! 

r. in 

= dv J d" v, S(i--J2 Vi) n exp (-v,) 



(13) 



(14) 



The predicted structure functions are 

= p 1 + 3?; 

^ 3 l+pT v 

which - although a nominally different form from Eq. [8] - gives 
nearly identical scalings for the choice T r w (7/6) | ln/3| (at least to 
P~20)Q 

4 For the "standard" SL assumptions in incompressible turbulence (/? = 
7 = 2/3), T v ft: 0.05 gives nearly identical structure functions. For the 



Boldyrev 1 20021 model often applied to super-sonic turbulence (j3 = 1/3, 
7 = 2/3), T r = 0.12 gives nearly identical structure functions. 
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In |Castaing| J1996) , the form above was motivated by ba- 
sic considerations of the most general form of an infinitely di- 
visible cascade. However, Yakhot (1998) independently obtained 
exactly the same scaling exponents (in the parameter B = l/T v 
derived therein) studying directly the symmetry properties of the 
Navier-Stokes equations. Moreover they derive T v w 1/20, rather 
than adopting it phenomenologically (for the case of nearly- 
incompressible sub-sonic turbulence, appropriate to match the orig- 
inal SL model; T r should be larger for compressible turbulence). 
They also directly derive the PDF of compressible (longitudinal) 
velocity fluctuations, which is exactly what we require to derive 
the PDF of density fluctuations. And indeed to first order in T v <C 1, 
the cascade model above captures the non-Gaussian behavior of the 
PDFs calculated in |Yakhot||l998| P1 

For our purposes here, this cascade model of the density PDF 
also has a plain interpretation in terms of thinking of the density and 
velocity fields as the product of a continuous-time multiplicative 
random relaxation process. As with the quantized case, the vari- 
able A represents some "number of events" over a dynamic range 
ln(?'2/ri) in scale. But we see that the function Gr in Eq. |13| is 
the convolution over a Poisson waiting time distribution for each 
quantized number of events. In other words, the PDF in Eq. [5]is a 
stationary result in a medium in which perturbations are driven by 
discrete, random events (point processes) with a constant average 
rate on a given scale, which then damp exponentially on a char- 
acteristic timescale (integrated over infinite timescales). T~ [ is a 
dimensionless "rate parameter"; for small values of T v ~ , events 
are widely spaced in time, so the variance between "levels of de- 
cay" at a given time is large, whereas for large values the variance 
is smaller (the medium is more smooth). The total variance of this 
distribution is 2XT V 2 , i.e. proportional to the "number of events" in 
some time, times the variance "per event." 

Chainais (20061 show that such a compound log-Poisson cas- 
cade corresponds to the most general form of an infinitely divisible, 
multidimensional cascade which can be decomposed into a combi- 
nation of point processes ("events") and random multipliers (for 
example, the Kolmogorov 1962 model, while ultimately producing 
some similar results, cannot be synthesized via a hierarchy of point 
processes except in one dimension). 

4.3 Testing this Explanation in Simulations 

We have already tested that the function in Eq. [5] provides an ex- 
ceptionally good fit to the simulation data. However we can further 
examine whether this might be connected with physical intermit- 
tency, as suggested by the model above. 

In Fig. [5] we show the velocity structure function scalings f„ 
predicted for pure self-similar turbulence, for the SL model with 
two commonly-adopted choices of 3 and 7 (Eq. 8). and for the 
Castaing ( 1996) model above with two values of 7^, (Eq. |14). We 
denote in this section the value T r as value in Eq.[l4]fit to the veloc- 
ity structure functions, to clearly discriminate it from the T in the 
fit to the density PDF. As noted above, for the appropriate choices 
of T v , the SL and Boldyrev (2002) structure functions are almost 
exactly reproduced, albeit with one fewer free parameter. We also 
compare the structure functions measured by |Pan et al.| ( |2009) from 
the simulations in Kritsuk et al. (2007 ), which agree very well with 
the |Boldyrev] ( |2002^ model as well as the T v = 0.12 |Castainglp996) 
model. Increasing T v represents systematically more intermittent 



turbulence, from T r — > (Kolmogorov) to T r — > oo (£„ — > 1 for all 
n > 1, as expected in pressure-free Burgers turbulence which is a 
random superposition of infinitely strong shocks). 

As noted above, the|Castaing (1996) r v -models give nearly 
identical structure functions to the SL model for the choice 



r v «-2 m /9 = --ln( l 1— 

6 6 V d-D a 



(15) 



5 To see this, insert the PDF of Eq.^nto Eq. 29 of |Yakhot|l998] dropping 
terms of 0{T, 2 ) ~ 1/400. 



where here d is the overall dimension of the system (d — 3 in every 
case here) and is the "fractal dimension" of the "most singular 
objects" (Doo — 1 for one-dimensional worms and lines, Doo = 
2 for sheet-like shock structures), while 7 is the singularity index 
(scaling of the size dependence of these structures; 7 = 2/3 for 
Kolmogorov, incompressible turbulence, while 7=1 for Burgers 
turbulence). 

We compile in Fig.[5]the velocity structure functions measured 
for all the simulations for which we have measured density PDFs, 
wherever available (see Table [TJ. For each, we fit a model of the 
form of Eq.[l4]to determine a best-fit T v (for convenience, assuming 
extended self-similarity). We see that there is a range of structure 
function shapes, but in every case a very good fit is possible with 
just one free parameter (T v ). 

Finally, we compare these T r values fitted to the velocity struc- 
ture functions, to the T determined from fitting the density PDF, for 
the same simulations. Unfortunately, the limited number of struc- 
ture function determinations means this plot has little data at large 
T. However there are a couple of interesting notes. At T > 0.1, 
we see T r ~ T, albeit with just a few data points and non-trivial 
scatter. At small T < 0.1, we see that essentially all the velocity 
structure functions asymptote to values consistent with the predic- 
tion for Boldyrev (2002|-like turbulence (T r ~ 0.12), even though 
the nominal T fit to the density distributions can be much smaller. 
But we caution here that at very small T the fractional uncertain- 
ties (fitting to the density PDF) are large (because the effects are 
very small), evident in the large error bars. Nevertheless there ap- 
pears to be a trend in the sense that the simulations which have the 
largest deviations from log-normality in the density PDF also ex- 
hibit the largest deviation from even Boldyrev (2002)-like scaling 
in the velocity structure functions - both with the sense of being 
"more intermittent." 

This is also consistent with the trend in Fig.|3]between T and 
S\ np or, similarly, compressive Mach number M c . 

5 DISCUSSION 

We have proposed an alternative (non-lognormal), physically moti- 
vated fitting function to describe the density distributions in ideal- 
ized simulations of isothermal turbulence. The proposed function, 
Eq.[5] provides a remarkably good fit to the data, over a huge range 
of Mach numbers and variance in the simulations, apparently in- 
dependent of the driving mechanism (solenoidal vs. compressive), 
presence or absence of magnetic fields, and numerical methods in 
the simulations. 

Although it is widely assumed that density PDFs in isother- 
mal turbulence are log-normal, we show that in fact a simple log- 
normal model is orders-of-magnitude less accurate as a description 
of these simulations, compared to the proposed function (especially 
in the wings of the distributions). Even higher order lognormal- 
skew-kurtosis models do little better. This manifests in a variety of 
ways, but particularly relevant is that the moments of the density 
distribution do not follow the predicted lognormal scalings l |3|4| ). 

This means, for example, that some previous studies which 
used the mean of In p as a proxy for the variance, or assumed that 
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the variance should be identical whether the PDF was volume- 
weighted or mass-weighted, could be significantly biased. The ro- 
bustness of the fits presented here allows us to recover the vari- 
ance 5 in the density distribution nearly identically whether we 
measure it directly from the numerical PDF, or treat it as a fit- 
ting parameter. This is itself an important check of numerical con- 
vergence in simulations (i.e. it suggests that there is not sufficient 
mass/volume in unresolved tails of the distribution to significantly 
bias the dispersion-Mach number relation and other correlations). 
But the variance fitted in a lognormal model often differs from that 
directly calculated from the simulation data, and can depend sig- 
nificantly on how the fit is weighted (e.g. whether equally in linear 
or log P(lnp)). 

Most simulations directly measure the variance used in the 
density dispersion-Mach number relation (without reference to a 
fitting function), but even in this case it is still important to spec- 
ify how the statistics are weighted; in Fig. [6] we show that quite 
different relations emerge depending on whether volume-weighted 
or mass-weighted dispersions are used in this correlation|^This dis- 
tinction is also critical to compare certain different numerical meth- 
ods (Lagrangian vs. Eulerian codes). If it were neglected, one might 
arrive at different values of, for example, the b parameter (relating 
A4 and M c ) in the density dispersion-Mach number relation de- 
pending on the density weighting and numerical method. 

Even given the same variance, it is clear that the proposed fit- 
ting function with non-zero T much more accurately describes the 
behavior of the "wings" of the density distribution. This may be 
especially important to models of star formation and fragmenta- 
tion in the ISM (see references in § [TJ, which depend on calcu- 
lating the mass or volume fraction above some threshold density 
(where e.g. self-gravity becomes important). If the threshold den- 
sity of interest is sufficiently far outside the "core" of the distri- 
bution, our comparisons here suggest it may be critical to account 
for the non-lognormal character of the density PDFs. In extreme 
cases, orders-of-magnitude errors could result from simply extrap- 
olating the density-Mach number relation and assuming a perfectly 
log-normal PDF (as is typically done in these models). 

We also provide a physical motivation for the proposed den- 
sity PDF. We discuss (§[T|l why mass conservation implies some 
asymmetry, breaking the conditions needed to drive the PDF to ex- 
actly lognormal; this asymmetry can be captured in the proposed 
PDF. And we specifically show that our proposed Eq.|5]can be de- 
rived from multi-fractal cascade models of the longitudinal (com- 
pressive) velocity increments, in a more general formulation of 
the well-known |She & Lev eque (1994) hierarchy owing to |Cas-| 
|taing| < | 1996| >. This suggests that the strength of deviations from 
lognormal may be directly linked to intermittency, describable by 
the same parameter as the intermittent velocity structure functions. 
And it does appear that the same deviations from self-similarity 



6 The mass-weighted variance S\ np ,M appears to scale as ln(l + as 
generally assumed. The volume- weighted variance S\ n „ y , however follows 
a steeper relation consistent with the trend we see in T (S) . Parameterized 
on its own, this appears to be closer to the dilogarithm Li^ (A4j5), where we 
can define 

r-i 



Stop,v«Ii2(M?)s / \n{\+M 2 c x) 
Jo 



dv 



(16) 



Interestingly, this is closer to the assumption that the total variance is the 
sum of the variance on all scales, where (scale-by-scale) the contribution to 
the variance is dSinp.vTO « ln(l + (M?[i?]))dlnR, as proposed in |Hop 
|kms||2012a| . 



appear directly in the structure functions of hip (see Kowal et al. 
|2007|[Liu~& Fang 2008). We also provide an alternative derivation 
based on a simple "counts in boxes" model in §[Aj this is heuristic 
(analogous to the discrete /? models used to describe some cases 
of the |She & Lev eque 1994) model, but it is instructive and lends 
a physical interpretation to the parameters in the fit, based on the 
degree of variance allowed "per step." In this scenario, the limit 
T — > oo corresponds to infinitely intermittent (Burgers) turbulence 
(a superposition of infinitely strong shocks), which would produce 
a highly non-normal density PDF, while T — > corresponds to a su- 
perposition of infinitely many smaller multiplicative perturbations, 
leading by the central limit theorem to a lognormal PDF. And in- 
deed we appear to see in the simulations that as the strength of the 
compressive component of the Mach number increases - i.e. as the 
fractional deviation induced by a strong shock, for example, be- 
comes larger - T increases systematically. In practice, the largest 
values of T we fit correspond to dimensions of intermittent struc- 
tures « 2.3, which agrees well with what directly measured in 
simulations of highly supersonic flows ( |Schmidt et al. { 2008 1. 

This provides an additional, powerful motivation to further 
study the deviations from log-normal behavior in density PDFs (in- 
dependent of the pure accuracy of the fitting function). The T pa- 
rameter here provides a direct diagnostic of the effects of intermit- 
tency on the density distribution in compressible turbulence, which 
remains poorly understood. To date, there has been very little ana- 
lytic work developing a framework in which to understand intermit- 
tency in the density field (as opposed to the velocity field), though 
simulations probing this have advanced tremendously. Measuring 
T alone, and understanding better how it scales with e.g. the Mach 
number of the simulations, can probe the nature of intermittent den- 
sity structures and how they depend on the compressibility of turbu- 
lence. Already here we see hints of a continuous change in the char- 
acter of these structures with J\A C - Comparing the T-values inferred 
from the density PDF to the full density and log-density structure 
functions provides a powerful and direct test of the cascade models 
described above. Related to this, understanding whether the density 
distribution does actually change shape with scale (in the manner 
predicted here), as is well-established for the velocity field, is fun- 
damental to any full understanding of the statistics of compressible 
turbulence. And comparison of T values and log-density structure 
functions to the velocity structure functions, expanding the lim- 
ited statistics available here, can test whether the same structures 
dominate the statistics of both the density and velocity fields. In 
particular, this is critical to understand whether something like a 
density-Mach number relation can be generalized scale-by-scale. 

That said, we strongly caution that some of the non-lognormal 
behavior in the PDFs fit here may be purely numerical. Because the 
high-density tail of the distribution is dominated by narrow, strong 
shocks and occupies little volume, the resolution demands are ex- 
treme, and convergence even at > 1000 3 cells/particles may not be 
complete (see Federrath et al. 2010 and references therein). Tech- 
nical details of the forcing scheme may bias the low-density tail 
as well (Konstandin et al.|2012) . And differences remain between 
different numerical methods (Price & Federrath 2010); it is interest- 
ing to note that in the one otherwise identical experiment we com- 
pile, SPH methods appear to show smaller T than fixed-grid meth- 
ods (i.e. are more biased to high versus low densities), which may 
fundamentally relate to the Lagrangian vs. Eulerian nature of the 
codes. So considerable caution is needed when giving any physical 
interpretation to the fitted T, especially for any single simulation. 
That said, significant deviations from lognormal behavior manifest 
here even where the PDFs appear to be numerically converged. And 
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in any case, as a pure phenomenological fitting function (indepen- 
dent of interpretation) the proposed function is a far more accurate 
description of the simulations than a lognormal. 

Finally, we also caution that in real astrophysical turbulence, 
much larger deviations from the ideal behavior discussed here 
should be expected for many reasons. Real systems are not exactly 
isothermal, and this leads to strongly non-normal density PDFs (see 
e.g. |Scalo et al.| 1 998 ; though in Appendix [B] we consider a gener- 
alized version of Eq. [3] that may be able to capture these features). 
Perhaps more important, self-gravity (e.g. Federrath & Klessen in 
preparation), large scale non-random perturbations in the ISM (spi- 
ral arms, mergers, etc.), and correlated feedback structures (radi- 
ation and SNe-induced "bubbles") can introduce very large devi- 
ations from Gaussianity (although in simulations, it seems that a 
roughly normal PDF may still apply to each such ISM phase sepa- 
rately; see e.g. Hopkins et al. 2012 1. In any case, we are not arguing 
that the corrections here are necessarily most important in the ISM, 
but they may be critical to compare models and simulations, and to 
understand the behavior of even ideal turbulence. 
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Table 1. Simulations Used in This Paper 



Reference 








Method 


Forcing 


Resolution 


M 


M c 


Ma 


S{np.V 


■^ln p.M 


T 


T r 


^rms 


u rms 


Federrath etal.{2010| 


UG 


C 


1024 3 


5.6 


5.5 


oo 


9.23 


3.23 


0.42 ±0.06 


0.54±0.10 


0.045 


0.77 


Federrath efaTl|2Ulo| 


UG 


s 


1024 3 


5.3 


1.8 


oo 


1.71 


1.45 


0.056 ±0.02 


0.17 ±0.03 


0.064 


0.42 


Schmidt et al. 1 2009) 


AMR 


c 


768 3 


2.9 


2.5 


oo 


3.33 


1.79 


0.23 ±0.03 


0.46 ±0.09 


0.055 


0.58 


Pnce"& Federrath |2010} 


SPH 


s 


512 3 


10 


3.3 


oo 


2.24 


2.01 


0.036± 0.014 


0.20 ±0.03 


0.14 


0.51 


Price & Federrath GOTO) 


UG 


s 


512 3 


10 


3.3 


oo 


2.97 


1.90 


0.16±0.02 


0.27 ± 0.07 


0.17 


1.23 


Kritsuk et al 


(20071 


UG/AMR 


M 


1024 3 


6 


2.2 


oo 


1.30 


1.28 


0.003 ±0.017 


0.12±0.01 


0.078 


0.11 


Konstandin et al. 


(20121 


UG 


c 


512 3 


15 


13.6 


oo 


12.93 


4.65 


0.38 ±0.1 6 


- 


0.057 


0.58 


Konstandin et al. 


poT2| 


UG 


s 


512 3 


15 


7.4 


oo 


4.58 


2.14 


0.25 ±0.08 


- 


0.11 


1.24 


Konstandin et al. 


(27JT2} 


UG 


c 


512 3 


5.5 


4.8 


oo 


8.17 


3.10 


0.34±0.16 


0.28 ±0.05 


0.091 


0.89 


Konstandin et al. 


pon) 


UG 


s 


512 3 


5.5 


2.4 


oo 


1.77 


1.47 


0.064 ±0.033 


0.16 ±0.04 


0.039 


0.60 


Konstandin et al. 


|2"0T2| 


UG 


c 


512 3 


2.0 


2.0 


oo 


5.51 


2.26 


0.29±0.12 


- 


0.049 


0.95 


Konstandin et al. 


pon) 


UG 


s 


512 3 


2.0 


0.95 


oo 


0.68 


0.62 


0.023 ±0.024 


- 


0.078 


0.36 


Konstandin et al. 


(27JT2") 


UG 


c 


512 3 


0.5 


0.4 


oo 


0.35 


0.27 


0.054 ±0.014 


- 


0.067 


0.74 


Konstandin et al. 


(Z0T2| 


UG 


s 


512 3 


0.5 


0.05 


oo 


0.019 


0.007 


0.018 ±0.003 


- 


0.21 


2.41 


Konstandin et al. 


|2"0T2| 


UG 


c 


512 3 


0.1 


0.1 


oo 


0.026 


0.023 


0.0001 ± 0.0002 


- 


0.12 


0.12 


Konstandin et al. 


pun) 


UG 


s 


512 3 


0.1 


0.0013 


oo 


7.4e-5 


6.9e-6 


0.002 ±0.001 


- 


0.43 


5.86 


Kowal et al. | 


2007 




UG 


s 


256 3 


7.1 


2.4 


7.1 


1.26 


1.23 


0.007 ± 0.045 


0.12 ±0.02 


0.050 


0.094 


Kowal et al. | 


2007 




UG 


s 


256 3 


2.34 


0.78 


7.4 


0.62 


0.55 


0.037 ±0.037 


0.10±0.03 


0.042 


0.21 


Kowafet al. | 


2007 




UG 


s 


256 3 


0.74 


0.25 


7.4 


0.064 


0.059 


0.023 ±0.009 


0.14 ±0.04 


0.14 


0.42 


Kowal"et al. | 


2007 




UG 


s 


256 3 


0.36 


0.12 


7.2 


0.0053 


0.0051 


0.015 ±0.003 


0.15 ±0.04 


0.24 


1.23 


Kowal"et al. | 


2007 




UG 


s 


256 3 


7.0 


2.0 


0.70 


1.39 


1.30 


0.012 ±0.045 


0.11 ±0.01 


0.081 


0.19 


Kowal et al. 


2007 




UG 


s 


256 3 


2.20 


0.63 


0.69 


0.66 


0.60 


0.034 ±0.037 


0.076 ±0.006 


0.10 


0.22 


Kowal et al. | 


2007 




UG 


s 


256 3 


0.68 


0.19 


0.65 


0.076 


0.068 


0.038 ±0.014 


0.074 ±0.006 


0.23 


0.66 


KowafetaTj 


2007 




UG 


s 


256 3 


0.33 


0.09 


0.64 


0.010 


0.0097 


0.024 ±0.005 


0.057 ±0.008 


0.29 


0.93 


Molina et al. 


120121 


UG 


M 


256 3 


17.6 


7.04 


oo 


3.68 


2.11 


0.32 ±0.07 




0.078 


0.57 


Molina et al. 


(2012 




UG 


M 


256 3 


16.8 


6.72 


7.0 


2.59 


1.90 


0.049 ±0.024 




0.13 


0.23 


Molina et al. 


pre 




UG 


M 


256 3 


10.6 


4.24 


oo 


2.89 


1.84 


0.27 ±0.08 




0.10 


0.58 


Molina et al. 


12012 




UG 


M 


256 3 


10.2 


4.08 


9.0 


2.16 


1.62 


0.049 ±0.030 




0.17 


0.29 


Molina et al. 






UG 


M 


256 3 


5.4 


2.16 


oo 


1.81 


1.37 


0.12 ±0.04 




0.075 


0.43 


Molina et al. 


hnn 




UG 


M 


256 3 


4.98 


1.99 


8.4 


1.35 


1.06 


0.076 ±0.026 




0.14 


0.43 


Molina et al. 


P0T2" 




UG 


M 


256 3 


2.21 


0.88 


oo 


0.59 


0.50 


0.058 ±0.022 




0.12 


0.48 


Molina et al. 


(20121 


UG 


M 


256 3 


2.09 


0.84 


8.1 


0.52 


0.43 


0.076 ±0.01 8 




0.17 


0.74 



(1) Reference: Source of simulation density PDF we fit (all data are taken directly from the published points in the PDFs). 

(2) Method: Numerical method. UG is uniform (Cartesian) grid; AMR is adaptive mesh refinement; SPH is smoothed-p article hydrodynamics . 

(3) Forcing: Turbulent velocity forcing scheme: C is purely compressive, S is solenoidal, M is mixed-mode driving (see |Federrath et al.|2010) . 

(4) Resolution: Numerical resolution. Refers to cell number for UG/AMR methods, particle number for SPH. 

(5) M: Mach number (at the driving scale), M = (vf) [ / 2 /c s . 

(6) M c '- Typical (rms) compressive Mach number, where quoted by the authors. Otherwise, we assume this = M for compressive forcing, = M /3 for solenoidal, and 
= 0.4.M for mixed-mode driving (see Federrath et al. 2010, Konstandin et al. 20121. With magnetic fields, the "effective" compressive component is smaller than this. 

(7) Ma- Alfen Mach number Ma = (vf) 1/1 /va- This is related to the ratio of thermal to magnetic pressure fi = P t h/-Pmag = 2c 2 /v A by /3 = 2 (Ma/M) 2 . 

(8) S]„ „ v : Volume-weighted variance in In (p) (computed directly from the simulation density distribution). 

(9) S{ np ,m- Mass-weighted variance in In (p) (directly from the simulation density distribution). This is equal to S\ np . v in the lognormal case, but smaller if T > 0. 

(10) T: Best-fit parameter T from Eq.|5] describing the deviations from a normal PDF. The 1 <r uncertainty is quoted, assuming all points in the PDF have equal 
logarithmic error bars (this is an upper lim it to the uncertainty fitted to any of the data sets shown). 

(11) T v : Best-fit parameter T v from Eq. | 14| describing the intermittency in the velocity structure functions, where available (note that the structure functions for the 
Federrath et al. 2010 simulations are actually in Schmidt et al. 20081. 

(12) 5^: Typical (rms) logarithmic deviation of the simulation data about the best-fit curve for the quoted T value (in dex). 

(13) (5^.: Deviations as Sj ms , but for the best-fit lognormal model (Eq. 2k. 
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APPENDIX A: A DERIVATION FROM A SIMPLE 
RANDOM COUNTS-IN-CELLS MODEL 

Consider the following simple, phenomenological model for 
stochastic variation in densities as a function of scale. Begin with a 
box containing total mass Mo and volume Vq, hence mean (volume- 
average) density po — Mo/Vo; without loss of generality choose 
units so Mq = Vo = po = 1. We then take a series of "steps." In 
each step, divide the volume V into two equal parts =V/2 (it does 
not matter how we perform this division). Select one of the vol- 
umes, and repeat until we reach a final volume Vf corresponding to 
the scale on which we wish to measure the densities. 

Consider two classes of steps. In a "normal" step, the density 
is smooth on the scale of interest, so each sub-volume contains ex- 
actly 1/2 the mass initially M within the volume V and density is 
conserved. In a "variable" step, there is some inhomogeneous struc- 
ture within the volume that makes the density non-uniform; con- 
sider the simple model where we assign to one sub-volume (V/2) 
a random fraction fi of the mass M (where ft is a uniform random 
variable between < fi < 1). By mass conservation, the other vol- 
ume has mass fraction /; = 1 — /;. But 1 — fi is also distributed as 
a uniform random variable between — 1. So for any final volume, 
selected at random, the mass is distributed as 2 - "' YYu^i fi> where 
m is the number of "normal" steps and n the number of "variable" 
steps. The final volume is = 2- (m+ " ) , so the final density is dis- 
tributed as p = 2" YYi=i fi (independent of m). So 

i=n 

ln(p)=nln2+^rn(/ i ) (Al) 

£=0 

i=n 

p ( u n = ^2HfO) =/ > (ln/,)®/ > (ln/ 2 )®...®P(ln/„) (A2) 

(A3) 



T(n) 



where ® denotes convolution. The latter equality follows because, 
for a uniform random variate / (P(f)df = d/, < / < 1), it is 
trivial to show P(lnf)d\nf = exp (-| ln/|) din/ for In/ < 0. So 
we just need the convolution over n independent exponentially dis- 
tributed variables, i.e. a gamma distribution. 

If the distribution of inhomogeneous structures is non- 
uniform, then the number n of "variable" steps is itself random; 
but assume they are statistically homogeneous on large scales, with 
a mean value (n) — A (for whatever final scale we average on) then 
this is just a counting process and n is Poisson-distributed. Com- 
bining, we obtain the final distribution 



\"e- x (-lnp + nln2)" 



(lnp-nln2) 



r(n) 



(A4) 



where the sum is restricted to lnp < n In 2 for each n. 

But this is clearly very closely related to our proposed Eq.|5] 
for the choice T — 1. In fact, Castaing ( 1996) showed that this is 
simply the "quantized" version of the PDF in Eq. [5] (i.e. has xa = 
ln2 > in the notation therein). This is because we have considered 
discrete steps in scale, with quantized "random" fluctuations (i.e. a 
binary choice between a multiplication by / or no change). If we 
instead allowed continuous steps in "scale" (i.e. differential steps 
in the volume element) but enforce the same "rate," per logarithmic 
change in volume, of encountering such a multiplicative "event" 
that divides the mass into / and 1 — fi, then we simply recover 



exactly the derivation in Castaing ( 1996) for T = 1 and, as a result, 

° AV- A (-lnp + A/2)'- 1 (lnp _ A/2) 

r(») 



(A5) 



for lnp < A/2 (/'(lnp) = otherwise). We can think of this as es- 
sentially taking the n "outside" of the sum (replacing it with A), 
since the number of "steps" in log-volume is always fixed (inde- 
pendent of the number of "random" steps) going from the box scale 
to the measurement scale. 

This is not a unique model for the density distribution, of 
course, but it by construction obeys mass conservation and exactly 
describes the density PDF on all scales, i.e. resolving the problems 
in §[T] We also stress that this is not the only model which arrives 
at the T = 1 density PDF; however it provides a useful context in 
which to interpret the PDF analogous to the use of the discrete fi 
model to motivate and interpret the intermittency model of |She &| 
|Leveque| ( [T994l >. 

We can further generalize this model to different T (at least 
for T <C 1). In a "variable" step, consider dividing the volume into 
T~ + 1 equal sub- volumes and taking all but one of them (a frac- 
tion T~ l /(T~ [ + 1) = 1/(1 + T)) as the "total sub- volume" in the 
step. Assume each sub-volume has a random fraction of the mass, 
subject to mass conservation as before, such that the mass-fraction 
PDF for any randomly-selected one of the sub-volumes is identical 
(again as before). For T <C 1, any such distribution leads to a total 
mass fraction in the T~ [ choisen sub- volumes distributed as f T to 
leading order (where / is again a uniform variate between — 1 ; 
this must be true because the PDF of the sum over 1 /T subvolumes 
- product of their individual probabilities - must converge to an 
/-like distribution). The final volume is then 2~ m (1 + T)~ n , and 
the final mass is distributed as 2 - "' II/ T , so the logarithm of this is 
distributed as 2"'" ~£ ln iff) = 2 ~"' T E m (fi)- Thus we obtain 



PQnp) 



r _1 [-lnp + ;iln(l + r)] 



(A6) 
(A7) 



which is again just the quantized version of our proposed density 
PDF (following Castaing 1996 to make this continuous with a con- 
stant "rate" parameter leads to exactly our proposed Eq.[5](. 

In this interpretation, then, A is equal to the mean number of 
"events" or "structures" producing density inhomogeneities - point 
processes in the cascade - integrated over a range of scales from 
the driving scale down to the averaging scale. To the extent that the 
density distribution converges, it is because these processes van- 
ish below some scale (becoming less significant below the sonic 
scale, and completely vanishing by the viscous scale). T represents 
the typical fractional variation associated with one such event (or 
equivalently the fraction of the mass on a given scale "associated 
with" a typical event). Because the density fluctuations associated 
with individual typical shocks, rarefactions, etc., are fractionally 
larger as the compressive Mach number increases, we naturally ex- 
pect a growing T with M c . 

APPENDIX B: A MORE GENERAL MODEL ALLOWING 
DIFFERENT SKEWNESS 

The cascade described above was derived by assuming a one-sided 
scaling in lnp. But as noted in Castaing (1996), mathematically 
speaking the same derivation could trivially generalize to fluctu- 
ations of opposite sign in lnp (though this may not be physical). 
Consider each fluctuation to have either positive or negative sign (in 
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log space) with probability of positive fluctuations p+. This gives 
a PDF: 

P{u)&u = y^ j P\{n) ^P + {m\n)P T {u\n,m)du (Bl) 

n— m— 

\ « — A 

A e 



Px(n) 



n 



P+(m n) = — -p+(l-p+) 

m ! (n — m ) ! 

\u\"~ 1 e~ M 

P T (u \n,m)= — — H (u \ n, m) 

(m— l)\(n — m— ij! 



H(u | n, m) 



f^dte-^'t"-"'" 1 (l+t)"" 1 (u > 0) 
J^dte-^'t'"-' (l+t)"-'"- 1 {u < 0) 



i r f +oc l 

«o(A, r) = In / du exp(T u)P(u)\ 

J— oo 

p = exp [r (;( + ;< )]po 

Here, Px is the Poisson probability of n events; P+ is the binomial 
probability of m "positive-biased" events (with per-event probabil- 
ity p+)\ Pt is the convolution of the exponential distribution of am- 
plitudes for a continuous-time relaxation process with m positive- 
amplitude events and n — m negative amplitude events; and m is 
determined by mass conservation: 



dupP(u)= duexp[T (u + u )]P(u) = 1 (B2) 



1 

Po . 

The difference between this and the formulation in the text is 
that here, the choice of p+ determines the level of skewness in the 
log-density PDF (p+ = corresponds to our model in the text). If 
we choose p+ — 1 we simply "mirror" the skewness in In p. The 
choice p+ = 1/2 is non-skew - although still non-Gaussian at a 
desired level set by T - since positive and negative perturbations 
are symmetric. This may be useful for fitting non-isothermal cases: 
for example, for a polytropic gas with c] oc p 7 ~' and 7 < 1, the 
PDF can be skewed towards higher densities (see |Scalo et al.|1998| 
Passot & Vazquez-Semadeni 1998 1, represented here by p+ > 1 /2. 
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